Napredne metode političkog istraživanja - Predavanje VIII - Korelacije i Linerana regresija
Fakultet političkih nauka, Podgorica
Olivera Komar, Nemanja Batrićević, Nemanja Stankov
April, 2021
Učitavanje već instaliranih paketa koje ćemo koristiti u analizi:
rm(list = ls()) # funkcija kojom čistimo radnu površinu
library('car') # Paket "car" ćemo koristiti za rekodiranje varijabli
library('psych') # Paket "psych" ćemo koristiti za prikazivanje mjera deskriptivne statistike
library('effects') # Paket "effects" ćemo koristiti za grafičku vizuelizaciju OLS regresija
library('dplyr') # Paket "dplyr" ćemo koristiti za sumiranje podataka
library('ggeffects') # Paket "ggeffects" ćemo koristiti za grafički prikaz interakcijskog efekta
library('ggplot2') # Paket "ggplot2" ćemo koristiti za grafički prikaz interakcijskog efekta
library('broom') # Paket "broom" ćemo koristiti za kreiranje tabele sa rezidualim, mjerenim i predviđenim vrijednostima regresijeUčitavanje dio datoteke Crnogorske izborne studije 2016.
mnes <- read.csv("mnes_2016_linreg.csv")
dim(mnes)[1] 1213 17head(mnes) lr god obr relig term_dps term_df term_ura term_dem efik intg nac pol
1 6 57 4 4 NA NA NA NA NA 5 CG Muski
2 NA 66 4 2 NA NA NA NA 4 5 CG Muski
3 NA 50 4 4 NA NA NA NA 5 4 CG Zenski
4 NA 70 4 3 NA NA NA NA 4 4 CG Muski
5 NA 66 2 2 NA NA NA NA NA NA CG Zenski
6 NA 68 2 3 NA NA NA NA NA NA CG Muski
polint zparl izl rur reg
1 4 Da Da Ruralno Sjeverna
2 2 Da Da Ruralno Sjeverna
3 4 Da Ne Ruralno Sjeverna
4 2 Da Da Ruralno Sjeverna
5 4 Da Ne Ruralno Sjeverna
6 4 Ne Ne Ruralno SjevernaBaza podataka sadrži informacije o ispitanicima i njihovim političkim preferencijama:
lr:pozicija ispitanika na skali lijevo-desnogod:godine ispitanikaobr:nivo obrazovanjarelig:stepen religionznostiterm_dps:partijski termometar: DPSterm_df:partijski termometar: DFterm_ura:partijski termometar: URAterm_dem:partijski termometar: DEMOKRATEefik:nema razlike u odnosu na to ko je na vlasti (spoljašnja efikasnost)intg:vjerovanje da se glasovi pošteno broje (izborni integritet)nac:nacionalna pripadnost ispitanikapol:pol ispitanikapolint:politička zainteresovanostzparl:stav prema što većoj zastupljenosti žena u parlamentuizl:izborna izlaznostrur:tip naselja u kojem ispitanik živireg:regija iz koje ispitanik dolazi
Korelacije
Možda je najčešće korišćeni pristup proučavanju odnosa između varijabli je korelacija. Postoje različite vrste korelacije, ali ovđe ćemo uglavnom govoriti o Pirsonovoj korelaciji, na koju se najčešće misli kada ljudi govore o korelaciji. Ona pokazuje povezanost između dvije intervalne varijable i implementira se u R kroz funkciju cor() i cor.test(). Prva jednostavno izračunava vrijednost koeficijenta korelacije, druga takođe izvodi statistički test kako bi vam rekao da li se korelacija statistički razlikuje od 0.
Funkcija cor() je korisna jer pruža mogućnost ispitivanja mnogih varijabli odjednom. Dakle, ako želimo ispitati sve korelacije između varijabli u datoteci koju koristimo:
cor(mnes[,1:10], use="pairwise.complete") lr god obr relig term_dps
lr 1.000000000 -0.03591071 -0.06605720 -0.06603875 0.03322555
god -0.035910707 1.00000000 -0.45656928 -0.11440266 -0.01655689
obr -0.066057202 -0.45656928 1.00000000 0.05140339 -0.01285286
relig -0.066038751 -0.11440266 0.05140339 1.00000000 0.01131959
term_dps 0.033225548 -0.01655689 -0.01285286 0.01131959 1.00000000
term_df 0.057693976 -0.04202890 -0.05013084 0.21605707 -0.26071008
term_ura -0.144728655 -0.16736524 0.07233529 0.08152652 -0.09488467
term_dem -0.089131609 -0.17633046 0.06940687 0.18089855 -0.19529052
efik 0.001983315 -0.05617688 0.06400845 0.10987867 0.09469162
intg -0.025496199 -0.03602413 0.05990342 -0.02817080 -0.70147492
term_df term_ura term_dem efik intg
lr 0.05769398 -0.14472866 -0.08913161 0.001983315 -0.02549620
god -0.04202890 -0.16736524 -0.17633046 -0.056176880 -0.03602413
obr -0.05013084 0.07233529 0.06940687 0.064008450 0.05990342
relig 0.21605707 0.08152652 0.18089855 0.109878667 -0.02817080
term_dps -0.26071008 -0.09488467 -0.19529052 0.094691616 -0.70147492
term_df 1.00000000 0.46599559 0.53063893 0.131271452 0.30917655
term_ura 0.46599559 1.00000000 0.57034656 0.183965590 0.18305092
term_dem 0.53063893 0.57034656 1.00000000 0.154684191 0.25060141
efik 0.13127145 0.18396559 0.15468419 1.000000000 -0.03439516
intg 0.30917655 0.18305092 0.25060141 -0.034395161 1.00000000Funkcija cor() je “probiljiva” kada su u pitanju nedostajuće vrijednosti i zato moramo naglasiti da želimo ispustiti sve slučajeve s vrijednostima koje nedostaju na varijablama za koju želimo izračunavanje korelacije. Argument "pairwise.complete" unutar funkcije zahtijeva od funkcije da koristi samo obzervacije koje imaju cjelovite podatke za obije varijable.
Viđeli smo da između velikog broja varijabli postoji korelacija. Vidimo da je korelacija između ideološke orijentacije i religioznosti veoma mala (-0.07). Možemo izvršiti i formalni test kako bi viđeli da li je ova korelacija statistički značajno drugačija od 0 ili nije. Nulta hipoteza u ovom slučaju kaže da ne postoji odnos između dvije varijable (korelacija je jednaka 0).
cor1 <- cor.test(mnes$lr, mnes$relig)
cor1
Pearson's product-moment correlation
data: mnes$lr and mnes$relig
t = -1.5351, df = 538, p-value = 0.1253
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.14958283 0.01844159
sample estimates:
cor
-0.06603875 Međutim, korelacija između ideloške pozicije i religioznosti nije daleko od statističke značajnosti, jer se interval pouzdanosti proteže preko 0.
Sa druge strane korelacija između starosne dobi i obrazovanja bi trebala biti značajna (-0.45) očito značajna.
cor2 <- cor.test(mnes$god, mnes$obr)
cor2
Pearson's product-moment correlation
data: mnes$god and mnes$obr
t = -17.799, df = 1203, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5001459 -0.4106861
sample estimates:
cor
-0.4565693 Korelacije se mogu grafički prikazati s dijagramom raspršenosti:
plot(mnes$god , mnes$obr)
Grafički prikaz nije u potpunosti jasan, iako se nazire šablon koji ukazuje da je korelacije negativna. Razlog zašto diagram raspršenosti ne izgleda toliko prijemčivo je taj što jedna od dvije varijable koje koristimo je zapravo ordinalna (kategorička) sa velikim brojem kategorija, zbog čega je tretiramo kao intervalnu. Ipak, to znači da će se veliki broj obzervacija (tačkica na grafiku) preklapati. Postoji više načina da se ovo ispravi i da dobijemo jasniju sliku. Kako nas trenutno zanima samo da prepoznamo korelaciju vizuelno, najednostavniji način je da izvučemo prosječne vrijednosti zavisne varijable (obrazovanje) za svaku o kategorija nezavisne (godine). Ovo možemo uraditi koristeći paket dplyr i funkciju summarise():
prosjeci <- group_by(mnes, god) #grupisanje podataka po starosnoj dobi
prosjeci_plot <- summarise(prosjeci, obr=mean(obr, na.rm=T)) # sumiranje prosječnog obrazovanja za svaku od prethodno grupisanih starosnih kategorija u novu tabelu
prosjeci_plot <- na.omit(prosjeci_plot) #isključivanje nedostajućih vrijednosti
prosjeci_plot # prikazivanje kreirane tabele# A tibble: 68 x 2
god obr
<int> <dbl>
1 18 4.83
2 19 4.12
3 20 4.89
4 21 4.62
5 22 5.58
6 23 4.9
7 24 6.06
8 25 5.4
9 26 5.61
10 27 6.05
# … with 58 more rowsplot(prosjeci_plot$god , prosjeci_plot$obr) #Grafičko prikazivanje sumiranih podataka
Kada grafički prikažemo prosječni nivo obrazovanja za svaku godinu starosti ispitanika u našem uzorku, odnos između dvije varijable se mnogo jasnije vidi.
Ukoliko želimo ovaj grafik uljepšati, možemo preimenovati nazive X i Y ose, dodati naslov, promijeniti boju i oblik tačaka…
plot(prosjeci_plot$god , prosjeci_plot$obr, # specifikovanje dvije varijable za koje smo računali korelacije
main = "Odnos između starosne dobi i obrazovanja", # izmjena naslova grafika
xlab = "Godine", # izmjena naziva X ose
ylab = "Nivo obrazovanja", # izmjena naziva Y ose
col = "navyblue", # izmjena boje tačaka
pch = 16, # izmjena oblika tačaka
cex = 1.5) # izmjena veličine tačaka
Napomena: Sirova vrijednost koeficijenta korelacije može zavarati i može ukazivati na jaču povezanost nego zapravo postoji. Za to bismo morali pogledati kvadrat koeficijenta korelacije, koji nam govori koliko se varijance u jednoj varijabli može objasniti varijansom u drugoj varijabli. Ako se koeficijent korelacije može kretati između -1 i 1, s 0 što ukazuje da nema povezanosti, korelacija od 0,5 možda se čini prilično snažnom povezanošću. Napokon, čini se da je to na pola puta između nikakve asocijacije i savršene asocijacije. Međutim, 0,5 samo znači da je objašnjeno 0,25 varijanse. Korelacija bi trebala biti oko 0,70 da bi mogli reći da je polovina varijanse na jednoj varijabli objašnjena vrijednostima na drugoj varijabli.
OLS Regresija
Možda najjednostavnija i najčešće korišćena analiza je linearna OLS (Ordinary Least Square) regresija. Postoje razne vrste regresije i osim OLS na ovom predmetu bavićemo se još i logističkog regresijom, ali u svakodnevnom govoru kada se kaže “regresija” misli se na OLS. Ona omogućava modelovanje intervalnih varijable kao linearne kombinacije (zbira) jedne ili nekoliko drugih intervalnih ili binarnih (dihotomnih) varijabli, tako da bismo na kraju imali grubu ideju o tome koliko bi se naša zavisna varijabla promijenila ako bi naša nezavisna varijabla promijenila za jednu jedinicu. Dakle, to je jednostavna, ali prilično fleksibilna i moćna statistička tehnika. Osnovni linearni model može se proširiti tako da obuhvati većinu analiza koje nam mogu pasti na pamet. Prosta regresija dobra je i zato što je relativno razumljiva. Njen osnovni princip je jednostavan - minimalizovanje zbroja kvadrata razlika između stvarnih i predviđenih vrijednosti - i može se prilično lako izračunati „ručno“.
Regresija je prikladna ako imamo intervalnu zavisnu varijablu, koja je više ili manje normalno distribuirana, intervalne ili binarne nezavisne varijable i razuman broj slučajeva koji su međusobno nezavisni. Odokativna pravila u odnosu na veličinu uzorka se razlikuju, ali sigurno ne bi bila dobra ideja raditi linearnu regresiju s manje od 30 slučajeva, posebno ukoliko ima mnogo nezavisnih varijabli. Regresija bi trebala biti “opravdana” i adekvatna ako ima više od 100 slučajeva i ne preveliki broj prediktora (nezavisnih varijabli). Što više želimo od podataka, tj. Što više stabilnih i procjena (estimates) želimo izvući iz podataka, to će nam više slučajve.
Koristeći datoteku s kojom smo se upoznali, pokušajmo modelovati ideološku poziciju koristeći varijable godine i obrazovanje.
U R-u linearni model može se kreirati s funkcijom lm(), koja ima iste argumente koje smo koristili u prethodnim funkcijama i analizama. Moramo odrediti formulu s zavisnom varijablom na lijevoj strani i nezavisnim varijablama na desnoj strani. U funkciji moramo odrediti i ime datoteke.
reg1 <- lm(lr ~ god + obr, data = mnes)
summary(reg1)
Call:
lm(formula = lr ~ god + obr, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-5.8733 -2.4752 -0.0355 2.7598 5.3249
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.871335 0.699632 9.821 <2e-16 ***
god -0.015404 0.009641 -1.598 0.1107
obr -0.172478 0.083622 -2.063 0.0396 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.246 on 573 degrees of freedom
(637 observations deleted due to missingness)
Multiple R-squared: 0.008694, Adjusted R-squared: 0.005234
F-statistic: 2.513 on 2 and 573 DF, p-value: 0.08195Prva stvar na koju bismo uvijek trebali obratiti pažnju je kvalitet modela (model fit). Pokazuju nam se dvije vrijednosti R-kvadrata na dnu izvoda iz modela. Od dvije vrijednosti uvijek bismo trebali gledati prilagođeni R-kvadrat, jer uzima u obzir i broj varijabli koje imamo u modelu i broj slučajeva kojima raspolažemo. Svaka dodatna varijabla u modelu, čak i ako uopšte nema asocijacije, povećava model fit na osnovu slučajnosti i ova mjera ispravlja tu procjenu.
Ovđe možemo viđeti da model ne odgovara podacimao, godine i obrazovanje pomažu da objasnimo otprilike 1% varijanse u ideloškom pozicioniranju. Takođe, vidimo da je 636 ispitanika izbrisano iz analize zbog nedostajućih vrijenosti. Veliki broj ispitanika je odbio ili nije razumio kako da se pozicionira na ideološkoj skali.
Probajmo u ovaj model dodati druge varijable. Recimo da očekujemo da će ispitanici koji su muškarci i koji žive na selu biti nešto više desno orijentisani od žena i onih koji žive u gradu. Varijable koje nijesu intervalne u linearnu regresiju uključujemo tako što od njih napravimo dihotomne varijable.
Prvo ćemo ih pretvoriti u vektore tipa karakter, kako bi ih lakše rekodirali.
mnes$rur <- as.character(mnes$rur)
mnes$pol <- as.character(mnes$pol)
mnes$nac <- as.character(mnes$nac)
mnes$reg <- as.character(mnes$reg)
mnes$izl <- as.character(mnes$izl)Prikaz kategorija po varijablama, kako bi jednostavnije mogli rekodirati varijable u sljedećem koraku:
table(mnes$rur)
Predgradje Ruralno Veci grad
503 438 271 table(mnes$nac)
AL BS CG SR
51 190 663 267 table(mnes$polint)
1 2 3 4
182 559 262 202 table(mnes$reg)
Centralna Juzna Sjeverna
498 274 441 table(mnes$izl)
Da Ne
1003 181 # Postoji više načina rekodirati varijablu. Vjerovatno najednostavniji način da se istovremeno rekodira veliki broj kategorija je korišćenjem funkcije recode() iz paketa "car". Međutim, često se dešava da postoje isti nazivi funkcija u više paketa. Ovo može kreirati probleme za samu analizu. U ovom slučaju postoji funkcija recode() i u paketu "dplyr" koji smo ranije koristili kako bi kreirali intuitivniji grafik za korelacije.
detach("package:dplyr", unload=TRUE) # Uklanjanje paketa "dplyr" nakon čega možemo bez problema rekodirati
mnes$grad_d <- recode(mnes$rur, "'Veci grad'=1; else=0")
mnes$zenski <- recode(mnes$pol, "'Zenski'=1; 'Muski'=0")
mnes$polint_d <- recode(mnes$polint, "1:2=0;3:4=1")
mnes$cg <- recode(mnes$nac, "'CG'=1;else=0")
mnes$sr <- recode(mnes$nac, "'SR'=1;else=0")
mnes$sjever <- recode(mnes$reg, "'Sjeverna'=1;else=0")
mnes$izl_d <- recode(mnes$izl, "'Da'=1;else=0")Pogledajmo da li je rekodiranje bilo uspješno?
head(mnes) lr god obr relig term_dps term_df term_ura term_dem efik intg nac pol
1 6 57 4 4 NA NA NA NA NA 5 CG Muski
2 NA 66 4 2 NA NA NA NA 4 5 CG Muski
3 NA 50 4 4 NA NA NA NA 5 4 CG Zenski
4 NA 70 4 3 NA NA NA NA 4 4 CG Muski
5 NA 66 2 2 NA NA NA NA NA NA CG Zenski
6 NA 68 2 3 NA NA NA NA NA NA CG Muski
polint zparl izl rur reg grad_d zenski polint_d cg sr sjever izl_d
1 4 Da Da Ruralno Sjeverna 0 0 1 1 0 1 1
2 2 Da Da Ruralno Sjeverna 0 0 0 1 0 1 1
3 4 Da Ne Ruralno Sjeverna 0 1 1 1 0 1 0
4 2 Da Da Ruralno Sjeverna 0 0 0 1 0 1 1
5 4 Da Ne Ruralno Sjeverna 0 1 1 1 0 1 0
6 4 Ne Ne Ruralno Sjeverna 0 0 1 1 0 1 0Sve izgleda u redu. Sada i dodatne varijable možemo uključiti u analizu:
reg2 <- lm(lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data= mnes)
summary(reg2)
Call:
lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg +
sr + sjever, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-7.1120 -2.1526 0.3648 2.3950 6.4494
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.903078 0.818657 10.875 < 2e-16 ***
god -0.018833 0.009615 -1.959 0.0506 .
obr -0.223437 0.091799 -2.434 0.0152 *
zenski -0.463734 0.273415 -1.696 0.0904 .
grad_d -2.363812 0.380120 -6.219 9.76e-10 ***
polint_d 0.092324 0.302920 0.305 0.7606
cg -0.520647 0.399384 -1.304 0.1929
sr -0.362316 0.406554 -0.891 0.3732
sjever -1.773211 0.364804 -4.861 1.52e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.124 on 565 degrees of freedom
(639 observations deleted due to missingness)
Multiple R-squared: 0.0942, Adjusted R-squared: 0.08138
F-statistic: 7.345 on 8 and 565 DF, p-value: 2.586e-09Ovaj model u cjelosti objašnjava 8% varijanse na zavisnoj varijabli. Dakle, 92% razlika između ispitanika u ideološkom pozicioniranju ostaje neobjašnjeno.
Presjek modela (intercept) nam govori da će kada ove varijable imaju vrijednost 0, onda će prosječna vrijednost zavisne varijable biti 8.9. Presjek nije uvijek instuitivan za interpretaciju, posebno ukoliko imamo veliki broj varijabli.
Kad je u pitanju koeficijen nagiba, interpretiramo isključivo efekte koji su statistički značajni. Takođe, veoma je važno voditi računa o tome kako je konkretna varijabla kodirana i koje kategorije služe kao referentne. Iz izvoda modela vidimo da statistički značajan efekat imaju varijable: obrazovanje, tip naselja i regija. Sa svakim novim stepenom obrazovanja ideološka pozicija pojedinca pomjeriće se u prosjeku za 0.22 poena ulijevo. Ovaj efekat je statistički značaj na nivou 0.05 (95% pouzdanosti). Takođe, ispitani koji žive u velikom gradu će u prosjeku biti pozicionirani za 2.36 poena lijevo u odnosu na one koji žive na selu ili u pregrađu gradova. Ovo možemo tvrditi sa 99.99% pouzdanosti da nijesmo napravili grešku I tipa. I na samom kraju osobe koje žive u sjevernoj regiji će u prosjeku biti za 1.77 poena pozicionirani u lijevo u odnosu na ispitanike iz južne ili centralne regije.
Pitanje:ako bi upoređivali efekat obrazovanja sa efektom regije, koji od dva ima snažniji uticaj na ideološko pozicionianje?
Napomena:interpretacija se radi uvijek uzimajući u obzir sve ostale varijable u modelu. Svaki pojedinačni koeficijen će vjerovatno biti drugačiji ukoliko bi promijenili model. Na primjer, probajmo iz modela isključiti varijablu regije. Što se mijenja?
reg2a <- lm(lr ~ god + obr + pol + grad_d + polint_d + cg + sr, data= mnes)
summary(reg2a)
Call:
lm(formula = lr ~ god + obr + pol + grad_d + polint_d + cg +
sr, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-6.6227 -2.3214 0.1167 2.4492 6.7023
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.01774 0.73521 9.545 < 2e-16 ***
god -0.02039 0.00980 -2.081 0.0379 *
obr -0.09402 0.08959 -1.049 0.2944
polZenski -0.44776 0.27881 -1.606 0.1088
grad_d -1.90538 0.37552 -5.074 5.29e-07 ***
polint_d 0.09012 0.30891 0.292 0.7706
cg 0.38876 0.35984 1.080 0.2804
sr -0.11334 0.41130 -0.276 0.7830
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.186 on 566 degrees of freedom
(639 observations deleted due to missingness)
Multiple R-squared: 0.05632, Adjusted R-squared: 0.04465
F-statistic: 4.826 on 7 and 566 DF, p-value: 2.698e-05Vidimo da je efekat tipa naselj značajno oslabio kada ne kontrolišemo za regiju iz koje dolazi ispitanik. Takođe, vidimo da sada starosna dob ispitanika ima značajan efekat, a da se efekat obrazovanja izgubio. Što to znači suštinski za našu teoriju?
Kao i u drugim prilikama, možda bi bilo korisnije i ovđe predstaviti rezultate grafički. Za regresije (i linearne modele općenito) to je olakšano paketom effects koji se može koristiti za izolovanje i crtanje efekta pojedine varijable zajedno sa intervalima pouzdanosti. Funkcija effect() izračunava marginalni efekt i koristi ime varijable i objekta modela, a generička funkcija plot() može se koristiti za grafičko prikazivanje efekta. Pogledajmo kako ovo izgleda za model koji smo upravo napravili.
plot(effect("obr", reg2),
main = "Uticaj obrazovanja na ideološku poziciju",
xlab = "Obrazovanje",
ylab = "Skala Lijevo-Desno")
Recimo da želimo testirati efekat istog modela na drugu zavisnu varijablu. Na primjer može nas zanimati kako nezavisne varijable utiču na percepciju političkih partija. U trećem regresionom modelu testiraćemo efekat ovih varijabli na to koliko se, na skali od 1 do 10, ispitanicima dopada DPS:
reg3 <- lm(term_dps ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data= mnes)
summary(reg3)
Call:
lm(formula = term_dps ~ god + obr + zenski + grad_d + polint_d +
cg + sr + sjever, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-6.8006 -2.6996 0.0933 3.3522 8.8958
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.116423 0.726711 9.793 < 2e-16 ***
god -0.004799 0.008315 -0.577 0.564
obr -0.095807 0.084896 -1.129 0.259
zenski 0.228806 0.244587 0.935 0.350
grad_d -0.201240 0.339187 -0.593 0.553
polint_d -1.550256 0.253672 -6.111 1.44e-09 ***
cg -0.416614 0.329539 -1.264 0.206
sr -4.384531 0.361081 -12.143 < 2e-16 ***
sjever 0.173170 0.331730 0.522 0.602
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.677 on 953 degrees of freedom
(251 observations deleted due to missingness)
Multiple R-squared: 0.1998, Adjusted R-squared: 0.1931
F-statistic: 29.75 on 8 and 953 DF, p-value: < 2.2e-16Dvije stvari su očigledno značajno različite. Prvo, vidimo da R-kvadrat koeficijent je značajno veći nego u prethodnom modelu (19%). Dakle, ovaj set varijabli mnogo bolje objašnjava stav prema ovo partiji, nego ideološku poziciju pojedinaca. Osim toga, vidimo da sasvim druge varijable imaju značajan efekat. Kako možemo interpretirati efekat nacionalne pripadnosti?
Objekat modela
Kad smo prethodno kreirali analize, sačuvali smo modele u objekte (poput reg1) i kada “pozovemo” ime objekta ili funkciju summary(), dobijamo neku vrstu izvoda ili pregleda onoga što se nalazi unutar objekta. Međutim, u stvarnosti objekat je složena lista koja sadrži što možemo zamisliti da nas zanima o modelu koji smo kreirali.
Pogledajmo objekt reg2 u kojem smo sačuvali rezultate prve regresije.
class(reg2)[1] "lm"Funkcija class() govori nam da imamo posla s objektom tipa lm - linearnim modelom. Ali, zapravo imamo posla s listom - onom istom koju smo ranije srijetali - ali s mnogo složenijom unutrašnjom strukturom. Unutrašnje djelove možemo viđeti ako na objektu koristimo funkciju str() koja će nam dati pregled strukture objekta.
str(reg2)List of 13
$ coefficients : Named num [1:9] 8.9031 -0.0188 -0.2234 -0.4637 -2.3638 ...
..- attr(*, "names")= chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ residuals : Named num [1:574] 1.27 5.54 0.34 -4.12 1 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ effects : Named num [1:574] -128.77 2.84 6.66 3.95 -15.78 ...
..- attr(*, "names")= chr [1:574] "(Intercept)" "god" "obr" "zenski" ...
$ rank : int 9
$ fitted.values: Named num [1:574] 4.73 4.46 5.66 6.12 7 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ assign : int [1:9] 0 1 2 3 4 5 6 7 8
$ qr :List of 5
..$ qr : num [1:574, 1:9] -23.9583 0.0417 0.0417 0.0417 0.0417 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:574] "1" "7" "9" "10" ...
.. .. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..- attr(*, "assign")= int [1:9] 0 1 2 3 4 5 6 7 8
..$ qraux: num [1:9] 1.04 1.05 1.07 1.06 1.02 ...
..$ pivot: int [1:9] 1 2 3 4 5 6 7 8 9
..$ tol : num 1e-07
..$ rank : int 9
..- attr(*, "class")= chr "qr"
$ df.residual : int 565
$ na.action : 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
$ xlevels : Named list()
$ call : language lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data = mnes)
$ terms :Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
$ model :'data.frame': 574 obs. of 9 variables:
..$ lr : int [1:574] 6 10 6 2 8 5 0 2 0 2 ...
..$ god : int [1:574] 57 63 25 48 26 54 19 54 21 37 ...
..$ obr : int [1:574] 4 5 8 4 4 4 7 4 7 7 ...
..$ zenski : num [1:574] 0 0 1 1 0 1 1 1 0 1 ...
..$ grad_d : num [1:574] 0 0 0 0 0 0 1 0 0 0 ...
..$ polint_d: num [1:574] 1 0 0 0 0 1 0 1 0 1 ...
..$ cg : num [1:574] 1 0 1 1 1 0 1 1 1 1 ...
..$ sr : num [1:574] 0 1 0 0 0 1 0 0 0 0 ...
..$ sjever : num [1:574] 1 1 0 0 0 0 0 0 0 0 ...
..- attr(*, "terms")=Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. .. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. .. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. .. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. .. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. .. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. .. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
..- attr(*, "na.action")= 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
.. ..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
- attr(*, "class")= chr "lm"Ovo izgleda prilično ružno i nečitko. Možemo probati da dobijemo manje informativni ali sažetiji pogled na imena različitih elemenata sadržanih u objektu. Možemo koristiti iste funkcije () koje smo ranije koristili za ispitivanje datoteka.
names(reg2) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "xlevels" "call" "terms"
[13] "model" Ovđe možemo pronaći sve vrste informacija o modelu, od koeficijenata modela do predviđenih vrijednosti zavisne varijable, grešaka i detalja o podacima i specifikaciji modela. Tako, na primjer, ako bismo htjeli izvući koeficijente iz modela, mogli bismo to učiniti upravo ovako:
reg2$coefficients(Intercept) god obr zenski grad_d polint_d
8.90307813 -0.01883318 -0.22343702 -0.46373399 -2.36381175 0.09232421
cg sr sjever
-0.52064731 -0.36231638 -1.77321086 Takođe, korisno je imati na umu da ako na objektu modelu korisite funkciju summary(), to stvara novu vrstu objekta,sa vlastitom unutrašnjom strukturom, koja bi za nas mogla imati korisne informacije.
summary_reg2 <- summary(reg2)
str(summary_reg2)List of 12
$ call : language lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data = mnes)
$ terms :Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
$ residuals : Named num [1:574] 1.27 5.54 0.34 -4.12 1 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ coefficients : num [1:9, 1:4] 8.9031 -0.0188 -0.2234 -0.4637 -2.3638 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
$ aliased : Named logi [1:9] FALSE FALSE FALSE FALSE FALSE FALSE ...
..- attr(*, "names")= chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ sigma : num 3.12
$ df : int [1:3] 9 565 9
$ r.squared : num 0.0942
$ adj.r.squared: num 0.0814
$ fstatistic : Named num [1:3] 7.34 8 565
..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
$ cov.unscaled : num [1:9, 1:9] 0.068661 -0.000557 -0.005408 -0.004924 -0.00336 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ na.action : 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
- attr(*, "class")= chr "summary.lm"Možemo viđeti da uglavnom sadrži iste informacije kao i objekt modela, poput koeficijenata i reziduala, ali i mnoge druge informacije. Pretpostavimo da smo kreirali više različitih modela modela i da želimo uporediti njihov kvalitet (model fit). To možemo jednostavno izdvojiti iz objekta sažetka.
summary_reg2$adj.r.squared[1] 0.08137558Iz objekata modela i izvoda (sažetka) tog objekta možete izvući puno više informacija. Imajte ovo na umu kada kreirate model i shvatite da su vam potrebne neke informacije o modelu kojih nema u sažetku ili ih ne možete tako lako izvući iz rezultata sažetka. Svakako sve informacije su objektu, samo ih nekad morate detaljnije potražiti.
Interakcija između nezavisnih varijabli
Posljednjim modelom pokušali smo procijeniti uticaj više nezavisnih varijabli na partijskii termometar. Cjelokupnim modelom mogli smo objasniti 19% varijanse u stavu prema DPS-u. Činjenica da model ne objašnjava najveći dio varijanse znači da bi se trebali zapitati da li je naš model dovoljno iscrpan? Da li postoje varijable koje imamo na raspolaganju koje birači uzimaju u obzir prilikom formiranja stava prema određenoj političkoj partiji?
U slučaju DPS, možemo očekivati da percepcija o tome koliko su izbori pošteni značajno utiče na stav prema ovoj partiji. Ukoliko je izborni integritet relevanatan faktor, možemo očekivati dvije stvari. Prvo, možemo očekivati da koeficijent bude statistički značajan; i drugo, da dođe do značajnog poboljšanja u kvalitetu modela (model fit - R kvadrat).
summary(reg3) # izvod iz prethodnog model
Call:
lm(formula = term_dps ~ god + obr + zenski + grad_d + polint_d +
cg + sr + sjever, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-6.8006 -2.6996 0.0933 3.3522 8.8958
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.116423 0.726711 9.793 < 2e-16 ***
god -0.004799 0.008315 -0.577 0.564
obr -0.095807 0.084896 -1.129 0.259
zenski 0.228806 0.244587 0.935 0.350
grad_d -0.201240 0.339187 -0.593 0.553
polint_d -1.550256 0.253672 -6.111 1.44e-09 ***
cg -0.416614 0.329539 -1.264 0.206
sr -4.384531 0.361081 -12.143 < 2e-16 ***
sjever 0.173170 0.331730 0.522 0.602
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.677 on 953 degrees of freedom
(251 observations deleted due to missingness)
Multiple R-squared: 0.1998, Adjusted R-squared: 0.1931
F-statistic: 29.75 on 8 and 953 DF, p-value: < 2.2e-16reg4 <- lm(term_dps ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever + intg, data= mnes) # Dopunjeni model
summary(reg4) # izvod dopunjenog modela
Call:
lm(formula = term_dps ~ god + obr + zenski + grad_d + polint_d +
cg + sr + sjever + intg, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-8.7591 -1.7382 0.2228 1.6108 8.8167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.345129 0.601368 17.203 < 2e-16 ***
god 0.003389 0.006704 0.505 0.613
obr -0.036742 0.067090 -0.548 0.584
zenski 0.317993 0.194112 1.638 0.102
grad_d 0.052991 0.266777 0.199 0.843
polint_d -0.846251 0.204275 -4.143 3.77e-05 ***
cg 0.214094 0.270397 0.792 0.429
sr -2.084357 0.303194 -6.875 1.19e-11 ***
sjever 0.000411 0.269931 0.002 0.999
intg -1.764975 0.069969 -25.225 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.787 on 868 degrees of freedom
(335 observations deleted due to missingness)
Multiple R-squared: 0.5424, Adjusted R-squared: 0.5377
F-statistic: 114.3 on 9 and 868 DF, p-value: < 2.2e-16Iz izvoda regresionog modela 4 možemo pročitati da negativan stav o izbornom integritetu vodi negativnoj percepciji DPS-a. Sa svakom jedinicom porasta nepovjerenja u izborni proces, stav prema DPS-u se smanjuje za 1.76 poena. Ovaj rezultat statistički je značajan na nivou p<0.001 (99.99%), što znači da imamo dovoljno dokaza da odbacimo nultu hipotezu koja glasi da ne postoji veza između ove dvije varijable. Osim toga, vidimo veoma značajan porast kvaliteta modela. Shodno prilagođenom R-kvadratu, dodavanjem samo jedne varijable procenat varijanse na zavisnoj varijabli koji možemo objasniti modelom porastao je sa 19 na čak 54.
Što se desilo sa ostalim varijablama? Možemo viđeti da su osim povjerenja u izborni integritet, kao i ranije, politička zainteresovanost i nacija relevantni faktori. Međutim, snaga efekta ovih varijabli sada je značajno manja nego ranije, otprilike prepolovljen. Što ovo znači? To znači da je makar dio varijanse zavisne varijable koji su ove nezavisne varijable objašnjavale zapravo dolazio od njihove povezanosti sa nepovjerenjem u izborni integritet. Iako je snaga efekta smanjena, sada smo sigurniji da je preostali uticaj na zavisnu varijabu zaista rezultat konkretnih varijabli a ne nekih nemjerenih fenomena.
Iako smo rekli da efekti pojedinačnih varijabli zavise od toga koje smo druge varijable uključili, modeli koje smo do ovog trenutka testirali pretpostavljaju da sve nezavisne varijable direktno utiču na nezavisnu. Međutim, istraživači često mogu pretpostaviti da određena varijabla može imati uticaj indirektno kroz drugu. Razlika između dva pblika modela teorijska, dok se statistički ovakavo odnos varijabli tretira kao interakcijski efekat. U R interakcija između dvije varijable specifikuje se koristeći * između dvije varijable.
Recimo da pretpostavljamo da postoji interakcije između dvije statistički neznačajne varijable: crnogorske nacionalnosti i nivoa obrazovanja. Mogli bi, na primjer pretpostaviti će se efekat nacionalnosti zavisiti od toga da li je osoba obrazovana ili ne.
reg4_int <- lm(term_dps ~ god + cg*obr + zenski + grad_d + polint_d + cg + sr + sjever + intg, data= mnes)
summary(reg4_int)
Call:
lm(formula = term_dps ~ god + cg * obr + zenski + grad_d + polint_d +
cg + sr + sjever + intg, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-8.6782 -1.8045 0.2315 1.6970 8.9455
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.044119 0.668588 16.519 < 2e-16 ***
god 0.003166 0.006687 0.473 0.6360
cg -0.972122 0.569268 -1.708 0.0881 .
obr -0.203044 0.097043 -2.092 0.0367 *
zenski 0.294374 0.193857 1.519 0.1293
grad_d 0.004895 0.266848 0.018 0.9854
polint_d -0.865933 0.203905 -4.247 2.40e-05 ***
sr -2.040114 0.302971 -6.734 3.01e-11 ***
sjever -0.088627 0.271837 -0.326 0.7445
intg -1.754754 0.069918 -25.097 < 2e-16 ***
cg:obr 0.267163 0.112912 2.366 0.0182 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.78 on 867 degrees of freedom
(335 observations deleted due to missingness)
Multiple R-squared: 0.5454, Adjusted R-squared: 0.5401
F-statistic: 104 on 10 and 867 DF, p-value: < 2.2e-16Na osnovu izvoda možemo viđeti da je zaista ovaj interakcijski efekat statistički značajan na nivou p < 0.05 (95% pouzdanosti). Efakt možemo interpretirati na sljedeći način: sa svakim dodatnim stepenom obrazovanja efekat pripadnosti crnogorskoj naciji (u odnosu na manjine) raste za 0.26 poena.
Ukoliko želimo efekat interacije grafički prikazati, možemo se poslužiti funkcijom ggpredict() iz paketa ggplot2:
grafik <- ggpredict(reg4_int, terms = c("cg", "obr"))
ggplot(grafik, aes(x, predicted, colour = group)) + geom_line()
Za slučaj da grafik nije previše čitljiv i da želite istaći samo neke od kategorija (linija) bilo iz estetskih ili suštinskih razloga, potrebno je samo u zagradu pored varijable dodati konkretne vrijednosti:
grafik_2 <- ggpredict(reg4_int, terms = c("cg", "obr[1,4,9]"))
ggplot(grafik_2, aes(x, predicted, colour = group)) + geom_line()
Na ovom grafiku vidimo tri regresione linije, pri čemu svaka označava zasebnu grupu Crnogoraca, shodno nivou obrzovanja. Možemo jasno viđeti da Crnogorci koji nemaju formalno obrazovanje imaju znatno više negativan u odnosu na pripadnike manjinskih naroda (crvena linija). Crnogorci koji imaju srednje obrazovanje imaju statistički jednaku percepciju DPS-a kao i manjinski narod (zelena linija). Na kraju, grupa visoko obrazovanih Crnogoraca ima značajno pozitivniji stav prema DPS-u u odnosu na manjinske narode (plava linija).
Provjera statističkih pretpostavki
Prije nego se pozabavimo regresionim pretpostavkama (dijagnostikom), poželjno je razjasniti dva ključna koncepta u regresionoj analizi: predviđene vrijednosti (fitted values) i rezidualne greške (resiual error). Oni su važni za razumijevanje dijagnostičkih grafika predstavljenih u nastavku.
Predviđene vrijednosti su vrijednosti zavisne varijable (Y) koje biste očekivali za zadate vrijednosti nezavisne varijable (X) prema izgrađenom modelu regresije (regresiona linija).
Iz dijagrama raspršenja vidi se da ne padaju sve tačke tačno na procijenjenu liniju regresije. To znači da se izmjerene vrijednosti mogu razlikovati od predviđenih vrijednosti. Razlika se naziva rezidualnim greškama, predstavljena je crvenim linijama.
U R možete lako dobiti podatke o predviđenim vrijednostima i rezidualima (ostacima) pomoću funkcije augment() iz paketa broom. Nazovimo izlazni reg.diag.metrics jer sadrži nekoliko važnih podataka korisnih za regresionu dijagnostiku.
reg.diag.metrics <- augment(reg4)U tabeli ima više varijabli:
Izmjerene vrijednosti svih varijabli .fitted: predviđene vrijednosti .resid: rezidualne greške
Pretpostavke
Linearna regresija pretpostavlja više stvari o podacima, kao što su:
- Linearnost. Pretpostavlja da je veza između prediktora (x) i zavisne varijable (y) linearna.
- Normalnost reziduala. Pretpostavlja da su rezidualne greško normalno distribuirane.
- Homogenost varijanse reziduala. Pretpostavlja da reziduali imaju konstantnu varijansu (homoscedasticity)
- Nezavisnost rezidualnih grešaka.
Ove pretpostavke ukoliko nijesu ispunjene mogu kreirati ozbiljne probleme za vjerodostojnost analize. Možemo ih provjeriti koristeći dijagnostičke grafike.
Dijagnostički grafici
Dijagnostički gafici se mogu dobiti koristeći osnovnu R funkcionalnost plot(), koja kreira 4 zasebna grafika koji odgovaraju ranije objašnjenim pretpostavkama. Ako u funkciju dodamo i redni broj onda ćemo dobiti svaki od grafika pojedinačno.
plot(reg4, 1) #Dijagnostički grafik 1
Grafik: Residuals vs Fitted.
Koristi se za provjeru pretpostavki linearnog odnosa. Vodoravna crta bez jasnih obrazaca pokazatelj je linearnog odnosa. U slučaju našeg modela, vidimo blaže odstupanje od ravne linije 0, na osnovu čega možemo zaključiti da je pretpostavka linearnosti nije narušena.
plot(reg4, 2) #Dijagnostički grafik 2
Grafik:Normal Q-Q.
Koristi se za ispitivanje da li su reziduali normalno distribuirani. Dobro je ako ostaci (reziduali) slijede ravnu isprekidanu liniju. Na ovom grafiku takođe vidimo blaža odstupanja na krajevima distribucije, ali ne bi se moglo reći da postoji značajno i sistematično odstupanje od regresione linije. Na osnovu ovoga možemo zaključiti da je su reziduali našeg modela normalno distribuirani.
plot(reg4_int, 3) #Dijagnostički grafik 3
Grafik:Scale-Location (or Spread-Location).
Koristi se za provjeru homogenosti varijanse reziduala (homoscedastičnost). Vodoravna linija sa jednako raširenim tačkama dobar je pokazatelj homogenosti.To nije skroz slučaj u našem primjeru, reklo bi se da imamo problem je varijansa veća u gornjim djelovima distribudije. Što može biti uzrok ovome?
plot(reg4, 4) #Dijagnostički grafik 4
Grafik:Residuals vs Leverage.
Koristi se za identifikaciju uticajnih slučajeva, odnosno ekstremnih vrijednosti koje mogu uticati na rezultate regresije kada se uključe ili isključe iz analize.
Četiri grafika prikazuju tri najekstremnije tačke podataka označene brojevima redova podataka u bazi podataka. Mogu biti potencijalno problematični. Možda biste ih trebali pažljivo pogledati kako biste provjerili postoji li nešto posebno za tu obzervaciju ili su to možda samo greške u unosu podataka.